function [posBinary, negBinary, Th, SessCorrData] = createBinaryMaskISCDifference(path, alpha)
% This function computes binary brain masks based on the thresholded ISC session comparison (ISC difference)
% maps.
%
% NOTE: This code computes binary masks only for ISC difference maps generated based on two sessions 
%       (i.e. one session comparison: Session 1 vs Session 2), and does not compute masks for time- 
%       and frequency-varying ISC maps.
%
% Usage: [posBinary, negBinary, Th, SessCorrData] = createBinaryMaskISCDifference(path, alpha)
%
% Inputs: path = string containing a full pathname to the setup-mat-file saved in the ISC analysis folder
%         alpha = a string containing a multiple comparison (FWER) corrected critical p-value. 
%                 Options are '0.05', '0.01' or '0.001'. If alpha is omitted, the default is '0.05'. 
%
% Outputs: posBinary = binary mask for the areas > Th when Sess1 > Sess2
%          negBinary = binary mask for the areas > Th when Sess2 > Sess1
%          Th = vector of thresholds for FWER alpha 0.05 0.01 and 0.001
%          SessCorrData = session corelation matrixes in (x,y,z,c) (Z values)
%
% Example: path = '/home/userX/ISCresults/Study1/StudyName.mat';
%          [posB, negB, Th, data] = createBinaryMaskISCDifference(path, '0.001')
%
% Juha Pajula, 17.9.2014

%check the number of inputs
if nargin==1
    alpha = '0.05';
end
%check the file existence
if exist(path,'file') == 0
   error('The path to setup file is not correct')
end
%check the threshold input
if strcmp(alpha,'0.05') 
    selTh = 1;
elseif strcmp(alpha,'0.01') 
    selTh = 2;
elseif strcmp(alpha,'0.001')
    selTh = 3;
else
    error('alpha must be string 0.05, 0.01 or 0.001!')
end

%load Params
load(path)

% load pointers:
load([Params.PublicParams.dataDestination 'memMaps'])

%aquire the session compare data (for the band 0):
SessCorrData=memMaps.PFMatMapSession.whole.band0.cor.sessComp1.Data.xyzc;
%compute the sum of Z statistics of Session correlation from the matrixes
avgSessCorrData = sum(SessCorrData,4);
negavgSessCorrData = (-1)*sum(SessCorrData,4);

%load thresholds
load([Params.PublicParams.dataDestination 'PFsession/ThPFBand0sessComp1win0.mat']) 

%positive binarymask from the results
posBinary = avgSessCorrData>Th(selTh);

%negative binarymask from the results
negBinary = negavgSessCorrData>Th(selTh);

%check the numbers:
disp(['Alpha: ' num2str(critVals)])
disp(['Th: ' num2str(Th)] )
disp('-------------------------------------------')
disp(['Used threshold with alpha: ' alpha ', value: ' num2str(Th(selTh))]);
disp(['Detected voxels when Sess1 > Sess2: ' num2str(sum(posBinary(:)))]);
disp(['Detected voxels when Sess2 > Sess1: ' num2str(sum(negBinary(:)))]);